#include "triangulator.h"
#include <iostream>
#include <cassert>
#include <algorithm>
using std::cerr;

#include <sstream>
string triangulator::trace_t(int t)
{
    if (t == -1)
        return "NULL";
    std::stringstream ss;
    ss << "T (" << pnts[tris[t].p[0]] << "), (" << pnts[tris[t].p[1]] << "), ("
        << pnts[tris[t].p[2]] << ")";
    return ss.str();
}

sign triangulator::is_inside(const triangle& t, const point& pp) const
{
    sign s1 = left_turn(pnts[t.p[0]], pnts[t.p[1]], pp);
    sign s2 = left_turn(pnts[t.p[1]], pnts[t.p[2]], pp);
    sign s3 = left_turn(pnts[t.p[2]], pnts[t.p[0]], pp);

    bool has_left = (s1 == POSITIVE || s2 == POSITIVE || s3 == POSITIVE);
    bool has_right = (s1 == NEGATIVE || s2 == NEGATIVE || s3 == NEGATIVE);

    return (has_left && has_right) ? OUT : INSIDE;
}

 sign triangulator::long_in_circle(const triangle& t, const point& pp) const
 {
     mpq_class ax = pnts[t.p[0]].x; mpq_class ay = pnts[t.p[0]].y;
     mpq_class bx = pnts[t.p[1]].x; mpq_class by = pnts[t.p[1]].y;
     mpq_class cx = pnts[t.p[2]].x; mpq_class cy = pnts[t.p[2]].y;
     mpq_class px = pp.x; mpq_class py = pp.y;

     mpq_class det1 = (ax - px) * (by - py) * ((cx * cx - px * px) + (cy * cy - py * py));
     mpq_class det2 = (ay - py) * (cx - px) * ((bx * bx - px * px) + (by * by - py * py));
     mpq_class det3 = (bx - px) * (cy - py) * ((ax * ax - px * px) + (ay * ay - py * py));

     mpq_class det4 = (cx - px) * (by - py) * ((ax * ax - px * px) + (ay * ay - py * py));
     mpq_class det5 = (bx - px) * (ay - py) * ((cx * cx - px * px) + (cy * cy - py * py));
     mpq_class det6 = (ax - px) * (cy - py) * ((bx * bx - px * px) + (by * by - py * py));

     mpq_class det = det1 + det2 + det3 - det4 - det5 - det6;

     if (det < 0)
         return OUT;
     if (det > 0)
         return INSIDE;
     return BOUND;
 }

 //copypast from Shewchuck
 sign triangulator::in_circle(const triangle& t, const point& d) const
 {
   double bdxcdy, cdxbdy, cdxady, adxcdy, adxbdy, bdxady;
   double alift, blift, clift;
   double det;
   double permanent, errbound;

   double ax = pnts[t.p[0]].x, ay = pnts[t.p[0]].y;
   double bx = pnts[t.p[1]].x, by = pnts[t.p[1]].y;
   double cx = pnts[t.p[2]].x, cy = pnts[t.p[2]].y;

   double adx = ax - d.x;
   double bdx = bx - d.x;
   double cdx = cx - d.x;
   double ady = ay - d.y;
   double bdy = by - d.y;
   double cdy = cy - d.y;

   bdxcdy = bdx * cdy;
   cdxbdy = cdx * bdy;
   alift = adx * adx + ady * ady;

   cdxady = cdx * ady;
   adxcdy = adx * cdy;
   blift = bdx * bdx + bdy * bdy;

   adxbdy = adx * bdy;
   bdxady = bdx * ady;
   clift = cdx * cdx + cdy * cdy;

   det = alift * (bdxcdy - cdxbdy)
       + blift * (cdxady - adxcdy)
       + clift * (adxbdy - bdxady);

   permanent = (abs(bdxcdy) + abs(cdxbdy)) * alift
             + (abs(cdxady) + abs(adxcdy)) * blift
             + (abs(adxbdy) + abs(bdxady)) * clift;
   double epsilon = std::numeric_limits<double>::epsilon();
   double iccerrboundA = (10.0 + 96.0 * epsilon) * epsilon;
   errbound = iccerrboundA * permanent;

   if (det < -errbound)
       return OUT;
   if (det > errbound)
       return INSIDE;
   else
       return long_in_circle(t, d);
 }

 template<typename T>
 vector<T>& operator +=(vector<T>& v1, const vector<T>& v2)
 {
     for (size_t i = 0; i < v2.size(); i++)
         v1.push_back(v2[i]);
     return v1;
 }

int triangulator::init(const std::vector<point> & points) 
{
    double max_x = 0;
    double max_y = 0;
    pnts = points;
    for (size_t i = 0; i < points.size(); i++)
    {
        if (std::fabs(points[i].x) > max_x)
            max_x = std::fabs(points[i].x);
        if (std::fabs(points[i].y) > max_y)
            max_y = std::fabs(points[i].y);
    }

    point p1 = point(-max_x * 100, max_y * 100);
    point p2 = point(0, -max_y * 100);
    point p3 = point(max_x * 100, max_y * 100);

    pnts += {p1, p2, p3};
    int tid = tris.size();
    tris.push_back(triangle(pnts.size() - 3, pnts.size() - 2, pnts.size() - 1, tid));
    in_triangulation.push_back(true);

    return tris.size() - 1;
}

int triangulator::localization(int t, int p)
{
    for (auto ct : tris[t].children)
        if (is_inside(tris[ct], pnts[p]) != OUT)
            return localization(ct, p);
//       WTF??

    return t;
}

int triangulator::reset(int t, int oldN, int newN)
{
    if (t == -1)
        return 0;
    for (int i = 0; i < 3; i++)
        if (tris[t].neigh[i] == oldN)
        {
            tris[t].neigh[i] = newN;
            return 1;
        }
    return -1;
}


void triangulator::create_new_triangles(int t, int p) 
{
    int p1 = tris[t].p[0];
    int p2 = tris[t].p[1];
    int p3 = tris[t].p[2];

    int t1i = tris.size();
    int t2i = tris.size() + 1;
    int t3i = tris.size() + 2;

    triangle t1(p1, p2, p, t1i);
    triangle t2(p2, p3, p, t2i);
    triangle t3(p3, p1, p, t3i);

    tris += {t1, t2, t3};
    in_triangulation += {true, true, true};

    tris[t1i].neigh = {tris[t].neigh[0], t2i, t3i};
    tris[t2i].neigh = {tris[t].neigh[1], t3i, t1i};
    tris[t3i].neigh = {tris[t].neigh[2], t1i, t2i};

    in_triangulation[t] = false;

    assert (reset(tris[t].neigh[0], t, t1i) != -1);
    assert (reset(tris[t].neigh[1], t, t2i) != -1);
    assert (reset(tris[t].neigh[2], t, t3i) != -1);

    assert (in_triangulation.size() == tris.size());

    tris[t].children += {t1i, t2i, t3i};

    check_triangle_and_flip(t1i);
    check_triangle_and_flip(t2i);
    check_triangle_and_flip(t3i);
}

void triangulator::check_triangle_and_flip(int t)
{
    if (!in_triangulation[t])
        return;
    for (int i = 0; i < 3; i++)
        flip(t, tris[t].neigh[i]);
}

void triangulator::flip(int t1, int t2) 
{
    if (t1 == -1 || t2 == -1)
        return;
    if (!in_triangulation[t1] || !in_triangulation[t2])
        return;

    assert (tris[t1].neigh[0] == t2 || tris[t1].neigh[1] == t2 || tris[t1].neigh[2] == t2);
    assert (tris[t2].neigh[0] == t1 || tris[t2].neigh[1] == t1 || tris[t2].neigh[2] == t1);

    auto in = [&](int p, int t) {
        return tris[t].p[0] == p || tris[t].p[1] == p || tris[t].p[2] == p;
    };

    int i1;
    for (i1 = 0; i1 < 3; i1++)
    {
        int nxt = (i1 + 1) % 3;
        if (in(tris[t1].p[i1], t2) && in(tris[t1].p[nxt], t2))
            break;
    }

    int p1 = tris[t1].p[i1];
    int p2 = tris[t1].p[(i1 + 1) % 3];
    int p3 = tris[t1].p[(i1 + 2) % 3];

    int i2;
    for (i2 = 0; i2 < 3; i2++)
        if (tris[t2].p[i2] != p1 && tris[t2].p[i2] != p2)
        {
            i2 = (i2 + 1) % 3;
            break;
        }
    int p4 = tris[t2].p[(i2 + 2) % 3];

    bool need_to_flip = false;

    sign s1 = in_circle(tris[t1], pnts[p4]);
    sign s2 = in_circle(tris[t2], pnts[p3]);

    if (s1 == INSIDE || s2 == INSIDE)
        need_to_flip = true;
    if (!need_to_flip)
        return;
    in_triangulation[t1] = false;
    in_triangulation[t2] = false;

    int nt1i = tris.size();
    int nt2i = tris.size() + 1;

    triangle nt1(p3, p4, p2, nt1i);
    triangle nt2(p4, p3, p1, nt2i);

    tris += {nt1, nt2};
    in_triangulation += {true, true};

    tris[nt1i].neigh = {nt2i, tris[t2].neigh[(i2 + 2) % 3], tris[t1].neigh[(i1 + 1) % 3]};
    reset(tris[t2].neigh[(i2 + 2) % 3], t2, nt1i);
    reset(tris[t1].neigh[(i1 + 1) % 3], t1, nt1i);

    tris[nt2i].neigh = {nt1i, tris[t1].neigh[(i1 + 2) % 3], tris[t2].neigh[(i2 + 1) % 3]};
    reset(tris[t1].neigh[(i1 + 2) % 3], t1, nt2i);
    reset(tris[t2].neigh[(i2 + 1) % 3], t2, nt2i);

    tris[t1].children += {nt1i, nt2i};
    tris[t2].children += {nt1i, nt2i};

    in_triangulation[t1] = false;
    in_triangulation[t2] = false;

    check_triangle_and_flip(nt1i);
    check_triangle_and_flip(nt2i);
}

    void triangulator::get_result(std::vector<triangle>& result, int t)
    {
        if (in_triangulation[t])
            result.push_back(tris[t]);
        visited[t] = true;
        for (int ct : tris[t].children)
            if ((ct != -1) && (!visited[ct]))
                get_result(result, ct);
    }

    std::vector<triangle> triangulator::delaunay_triangulation(std::vector<point> points)
    {
        random_shuffle(points.begin(), points.end());
        int init_triangle = init(points);
        in_triangulation[init_triangle] = false;

        for (int id = 0; id < (int)points.size(); id++)
        {
            int t = localization(init_triangle, id);
            create_new_triangles(t, id);
        }

        std::vector<triangle> result;
        visited.resize(tris.size(), false);
        get_result(result, init_triangle);
        std::vector<triangle> filtered;
        for (auto t : result)
        {
            if (t.p[0] >= (int)points.size() || t.p[1] >= (int)points.size() || t.p[2] >= (int)points.size())
                continue;
            filtered.push_back(t);
        }
        return filtered;
    }

point triangulator::circumscribed_circle(int t) const
{
    point p1 = pnts[tris[t].p[0]];
    point p2 = pnts[tris[t].p[1]];
    point p3 = pnts[tris[t].p[2]];

    double a1 = (p2.x - p1.x);
    double b1 = (p2.y - p1.y);
    double c1 = -0.5 * ((p2.x - p1.x) * (-p2.x - p1.x) + (p2.y - p1.y) * (-p2.y - p1.y));

    double a2 = (p3.x - p1.x);
    double b2 = (p3.y - p1.y);
    double c2 = -0.5 * ((p3.x - p1.x) * (-p3.x - p1.x) + (p3.y - p1.y) * (-p3.y - p1.y));

    double det = a1 * b2 - a2 * b1;
    double detx = c1 * b2 - c2 * b1;
    double dety = a1 * c2 - a2 * c1;

    return {detx / det, dety / det};
}

std::vector<std::pair<point, point> > triangulator::build_diagram(const vector<point>& points)
{
    std::vector<triangle> delaunay = delaunay_triangulation(points);
    visited = vector<bool>(tris.size(), false);

    std::vector<std::pair<point, point> > result;
    for (auto t : delaunay)
    {
        for (auto ti : t.neigh)
        {
            if (ti == -1)
                continue;
            if (!visited[ti])
                result.push_back(std::make_pair(circumscribed_circle(t.id), circumscribed_circle(ti)));
        }
        //assert (t.id < visited.size());
        visited[t.id] = true;
    }

    return result;
}
